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ESTIMATION OF CONSTANT AND TIME- VARYING DYNAMIC 
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Modeling viral dynamics in HIV/ AIDS studies has resulted in a 
deep understanding of pathogenesis of HIV infection from which novel 
antiviral treatment guidance and strategies have been derived. Viral 
dynamics models based on nonlinear differential equations have been 
proposed and well developed over the past few decades. However, it 
is quite challenging to use experimental or clinical data to estimate 
the unknown parameters (both constant and time-varying parame- 
ters) in complex nonlinear differential equation models. Therefore, 
investigators usually fix some parameter values, from the literature 
or by experience, to obtain only parameter estimates of interest from 
clinical or experimental data. However, when such prior information 
is not available, it is desirable to determine all the parameter esti- 
mates from data. In this paper we intend to combine the newly de- 
veloped approaches, a multi-stage smoothing-based (MSSB) method 
and the spline-enhanced nonlinear least squares (SNLS) approach, to 
estimate all HIV viral dynamic parameters in a nonlinear differen- 
tial equation model. In particular, to the best of our knowledge, this 
is the first attempt to propose a comparatively thorough procedure, 
accounting for both efficiency and accuracy, to rigorously estimate 
all key kinetic parameters in a nonlinear differential equation model 
of HIV dynamics from clinical data. These parameters include the 
proliferation rate and death rate of uninfected HIV-targeted cells, 
the average number of virions produced by an infected cell, and the 
infection rate which is related to the antiviral treatment effect and 
is time-varying. To validate the estimation methods, we verified the 
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identifiability of the HIV viral dynamic model and performed simu- 
lation studies. We applied the proposed techniques to estimate the 
key HIV viral dynamic parameters for two individual AIDS patients 
treated with antiretroviral therapies. We demonstrate that HIV vi- 
ral dynamics can be well characterized and quantified for individual 
patients. As a result, personalized treatment decision based on viral 
dynamic models is possible. 

1. Introduction. The AIDS pandemic continues to pose a serious threat 
to public health worldwide. AIDS has killed an estimated 2.9 million peo- 
ple, with close to 4.3 million people newly infected with HIV in 2006 alone. 
The total number of people living with HIV has reached its highest level 
(WHO web site). Although highly active antiretroviral therapy (HAART) 
regimens are effective in suppressing plasma HIV RNA levels (viral load) 
below the limit of detection, many patients may fail HAART treatments 
due to drug resistance, lack of potency, poor drug adherence, pharmacoki- 
netic problems and adverse effects. In addition, the complexity of regimens 
and lack of full understanding of the pathogenesis of HIV infection also 
pose great challenges to AIDS researchers. Over the past 2 decades, many 
mathematicians and statisticians have developed mechanism-based models 
and statistical approaches to assist in understanding HIV pathogenesis and 
have made significant contributions in this area [Ho et al. (1995), Wei et al. 
(1995), Perelson et al. (1996, 1997), Wu et al. (1999)]. In particular, differ- 
ential equation models have been widely used in describing dynamics and 
interactions of HIV and the immune system. Some survey on these models 
can be found in Perelson and Nelson (1999), Nowak and May (2000), and 
Tan and Wu (2005). 

For the mechanism-based models of HIV infection, one critical question 
is how to use experimental or clinical data to estimate the parameters in 
the nonlinear differential equation models which do not have closed-form 
solutions in most cases. Researchers have made a substantial effort to get an 
approximate closed- form solution under various assumptions, and then use 
the standard regression approach to estimate the dynamic parameters in the 
models [Ho et al. (1995), Wei et al. (1995), Perelson et al. (1996, 1997), Wu 
et al. (1999), Wu and Ding (1999), Wu (2005)]. But these approximations 
and assumptions may not always hold, in particular, for patients undergoing 
long-term treatment. 

In this paper we consider the following HIV dynamic model for patients 
under long term treatments [e.g., Perelson and Nelson (1999)]: 

(1.1) jT v (t) = A - P Tu(t) - V (t)Tu(t)V(t), 

(1.2) ^-Tj(t) = V (t)Tu(t)V(t) - 5Tj(t), 
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(1.3) -V(t) = N5Tj(t)-cV(t), 

at 

where Tjj is the concentration of uninfected target CD4+ T cells, Tj the 
concentration of infected cells, V the viral load, A the proliferation rate of 
uninfected target cells, p the death rate of uninfected target cells, 77 (t) the 
time-varying infection rate depending on antiviral drug efficacy, 5 the death 
rate of infected cells, c the clearance rate of free virions, and N the number 
of virions produced by a single infected cell on average. In this system, Tjj(t), 
Ti(t) and V(t) are state variables, and (\,p,N,5,c,r](t)) T are unknown dy- 
namic parameters. Notice that here we do not distinguish the infected cells 
by various subpopulations such as productively infected cells, latently in- 
fected cells and long-lived infected cells [Perelson et al. (1996, 1997)], since 
we intend to model the viral dynamics for an HIV-infected patient under 
long-term treatment. We also use a time-varying parameter rj(t) to model 
the infection rate since the infection rate may change nonparametrically due 
to the variation in treatment effect over time. Generally speaking, all kinetic 
parameters in this model could be time- varying [not only rj(t)]; however, it 
is usually unfeasible to do so due to the limited data collected in the clinical 
studies. Thus, this model provides a flexible although simple approach for 
studying long-term viral dynamics. In clinical trials or clinical practice, viral 
load, V(t) and total CD4+ T cell count, T{t) = T v {t) +Tj(t), are closely 
monitored and measured over time. 

In practice, some parameters can be fixed from previous studies and only 
the remaining parameters are needed to be estimated. However, when such 
prior knowledge is not available, it is important to estimate all viral dy- 
namic parameters, (X,p,N,5,c,r](t)) T , for each individual patient from the 
clinical measurements of V(t) and T(t), since the estimated dynamic param- 
eters may be used to guide clinical decisions for individualized treatment. 
Although AIDS investigators have tried to estimate some of these dynamic 
parameters based on viral load data [Ho et al. (1995), Wei et al. (1995), Perel- 
son et al. (1996, 1997), Wu et al. (1999)], none of them have successfully 
estimated all these dynamic parameters directly for an individual patient. 
Huang and Wu (2006) and Huang, Liu and Wu (2006) have made an effort 
to use the Bayesian method to estimate all these parameters, but strong 
priors are required for most parameters in order to make all parameters 
identifiable. In the work of Xia (2003), Filter, Xia and Gray (2005), Gray 
et al. (2005), and Ouattara, Mhawej and Moog (2008), all model parame- 
ters were estimated for HIV dynamic models with only constant coefficients. 
In this paper, to the best of our knowledge, this effort represents the first 
attempt to simultaneously estimate all the constant and time-varying HIV 
viral dynamic parameters for individual patients using both viral load and 
total CD4+ T cell count data. 
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In this paper we propose two estimation approaches. The first one is a 
multistage smoothing-based (MSSB) approach. We derive the direct rela- 
tionships between the unknown dynamic parameters and measurement vari- 
ables from the original differential equation models, and then we employ 
these relationships to formulate regression models for the unknown parame- 
ters. These regression models involve not only the time function of measure- 
ment variables, but also their derivatives. We propose using a nonparametric 
smoothing method, such as the local polynomial smoothing, to obtain the 
measurement variables and their derivatives which are substituted into the 
regression models to estimate the unknown dynamic parameters (includ- 
ing the time- varying parameter) in the second step. This approach avoids 
directly solving the ODEs by numerical methods and is computationally 
efficient. The second approach that we propose is to use a spline-enhanced 
nonlinear least squares (SNLS) method. This approach is to use splines to 
approximate the time-varying parameter so that the original ODE model 
becomes a model only containing constant parameters. The standard NLS 
method can be employed to estimate the unknown dynamic parameters and 
spline coefficients. This approach needs to use a numerical method to repeat- 
edly solve the nonlinear ODEs for the high-dimensional NLS optimization 
problem which is computationally challenging. Some cutting-edge computa- 
tional techniques are necessary to solve this problem. Note that the MSSB 
approach is computationally efficient, but the derived estimates are quite 
rough; however, these estimates can be used as the initial values for the 
SNLS method so that the two approaches can be efficiently combined for 
practical applications. 

The rest of this paper is organized as follows. In Section 2 we briefly 
discuss the identifiability of the HIV dynamic model, then propose the mul- 
tistage smoothing-based (MSSB) approach. In Section 3 we introduce the 
spline-enhanced nonlinear least squares (SNLS) method. In particular, a 
hybrid optimization technique combining a gradient method and a global 
optimization algorithm will be used to tackle the high-dimensional NLS op- 
timization problem. The simulation studies are presented to evaluate the 
performance of the proposed approaches in Section 4. We apply the pro- 
posed methods to estimate HIV viral dynamic parameters (including the 
time-varying infection rate) from data for two AIDS patients in Section 5. 
We demonstrate that some of these dynamic parameters are estimated from 
clinical data for the first time. We conclude the paper with some discussions 
in Section 6. The details of the proposed hybrid optimization algorithm are 
given in the Appendix. 

2. A multistage smoothing-based (MSSB) approach. Before introduc- 
ing the estimation methods, it is critical to determine whether all unknown 
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model parameters can be uniquely and reliably identified from a given sys- 
tem input and measurable output. The technique to answer this question is 
called identifiability analysis. There are mainly two types of identifiability 
problems: structural (mathematical) identifiability and practical (statistical) 
identifiability. Structural identifiability analysis techniques use model struc- 
ture information only (without knowledge of experimental observations) to 
determine whether all unknown parameters are uniquely identifiable. There- 
fore, it is also called the prior identifiability analysis. In contrast, practical 
identifiability analysis techniques should be applied after model fitting to 
verify the reliability of estimates, so it is called the posterior identifiability 
analysis [Rodriguez- Fernandez, Egea and Banga (2006)]. A thorough dis- 
cussion and review of identifiability techniques and theories is out of the 
scope of this study; for more details, the interested reader is referred to 
Ritt (1950), Bellman and Astrdm (1970), Kolchin (1973), Pohjanpalo (1978), 
Cobelli, Lepschy and Jacur (1979), Walter (1987), Vajda, Godfrey and Ra- 
bitz (1989), Ollivier (1990), Chappel and Godfrey (1992), Ljung and Glad 
(1994), Audoly et al. (2001), Xia and Moog (2003), Jeffrey and Xia (2005), 
Wu et al. (2008) and Miao et al. (2008). In this study we employed the 
differential algebra approach [Ljung and Glad (1994)] and it can be easily 
verified that all constant and time- varying parameters in models (1.1)-(1.3) 
are structurally identifiable. 

Parameter estimation for ODE models has been investigated using the 
least squares principle by mathematicians [Hemker (1972), Bard (1974), Li, 
Osborne and Prvan (2005)], computer scientists [Varah (1982)] and chemical 
engineers [Ogunnaike and Ray (1994), Poyton et al. (2006)]. Mathematicians 
have focused on the development of efficient and stable algorithms to solve 
the least squares problem. Recently statisticians have started to develop var- 
ious statistical methods to estimate dynamic parameters in ODE models. 
For example, Putter et al. (2002), Huang and Wu (2006), and Huang, Liu 
and Wu (2006) have developed hierarchical Bayesian approaches to estimate 
dynamic parameters in HIV dynamic models for longitudinal data. Li et al. 
(2002) proposed a spline-based approach to estimate time-varying parame- 
ters in ODE models. Ramsay (1996) proposed a technique named principal 
differential analysis (PDA) for estimation of differential equation models [see 
a comprehensive survey in Ramsay and Silverman (2005)]. Recently Ram- 
say et al. (2007) applied a penalized spline method to estimate the constant 
dynamic parameters in ODE models. Chen and Wu (2008, 2009) and Liang 
and Wu (2008) proposed a two-step smoothing-based approach to estimate 
both constant and time-varying parameters in ODE models separately. 

In this section we adopt the ideas from Liang and Wu (2008) and Chen 
and Wu (2008, 2009) to use the smoothing-based approach to estimate all 
dynamic parameters, including both constant and time-varying parameters, 
in model (1.1)— (1.3) in three stages. Stage I is to smooth the noisy data 
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to estimate the state variables and their derivatives; Stage II is to estimate 
constant dynamic parameters (A, p, c); Stage III is proposed to estimate both 
constant dynamic parameters, 5 and N, and the time-varying parameter 
rj(t). We introduce these three stages in detail in the following subsections. 

2.1. Stage I: Local polynomial estimation of the state variables. In this 
subsection we briefly describe how we use the local polynomial regression 
technique to estimate the functions T(t) and V(t) and their derivatives. 
Consider a general situation that a state variable, X(t), is observed at 
n time points, (Yi, . . . , Y n ), that is, Yi = X(ti) + e% for % = 1, . . . , n, where 
(ei, C2, ■ ■ ■ , e n ) are independent measurement errors with mean zero. Assume 
that the fourth derivative of X(t) exists. For each given time point to, we 
approximate the function X(ti) locally by a pth-order polynomial, X(ti) ~ 
X(to) + (U - t )X^(t ) + • • ■ + xV)(t )(U - t y/pl ± E?=oCj-(<o)(*i " to) 3 , 
for ti, i = 1, . . . ,n, in a neighborhood of the point to, where Cj(to) — X^\to) 
for j = 0,1,..., p. Following the local polynomial fitting [Fan and Gijbels 
(1996)], the estimators jfM(t) of xM(t) (v = 0,1,2 in our case) can be 
obtained by minimizing the locally weighted least-squares criterion, 

n ( p \ 2 

i=l I j=o J 

where K{-) is a symmetric kernel function, Kh(-) = K(-/h)/h, and h is a 
proper bandwidth. Let Y = (Y\, . . . , Y n ) T denote T p> t as an n x design 
matrix whose ith row is (l,tj — t, . . ., (t{ — i) p ) T , and W[ an n x n diagonal 
matrix of kernel weights, that is, diag{A^(£i — t), . . . ,Kh(t n — t)}. Assum- 
ing that the matrix Tj jt W t Tp jt is not singular, a standard weighted least 

squares theory leads to the solution, C = (T^WtTp^^Tj^WtY, where 
p = 1,2,3. We use the local linear regression to estimate X(t), the local 
quadratic regression to estimate X'(t), and the local cubic regression to 
estimate X^ 2 \t). Consequently, using the above notation, the estimators 
X^ (t) can be expressed as 

X®(t) =£l +1) (Tj +q jW t T 1+q , t r 1 Tj +q jW t Y for g = 0,l,2, 

where £ q is the (q + 2) x 1 vector having 1 in the (q + 1) entry and zeros in 
the other entries. 

We apply the local polynomial smoothing technique to the viral load data, 
V(t) and total CD4+ T cell count data, T{t) = T v (t) + T/(i) to obtain the 
estimates of V(t), T(t), V'(t) and T'(t) that will be used in Stage II, and 
the estimate of (t) that will also be used in Stage III. 



PARAMETER ESTIMATION IN NONLINEAR ODE MODELS 



7 



2.2. Stage II: PsLS estimates. In this subsection we apply the approach 
proposed by Liang and Wu (2008) to estimate constant dynamic parameters 
(A,/3,c) in models (1.1)— (1.3). Notice that we can combine equations (1.1) 
and (1.2), and obtain 

j t [Tu(t)+Tj(t)} = A - P Tu(t) - STj(t). 

Recalling that T{t) = Tj(t) + T v (t) and substituting T v (t) = T(t) - T r (t) in 
the above equation, we obtain 

±T(t) = \-p[T(t)-T I (t)]-5T I (t). 
Prom the above equation, we obtain 

T I = ^- + -^—T+-^—T\ 
p — o p — o p — o 

where T' = dT(t)/dt. Substituting this expression into equation (1.3) and 
letting a = a x = ffi and a 2 = we have 

(2.1) V\t) = a + aiT(t) + a 2 T'{t) - cV(t), 

where V(t) and T(t) for t = t\, t 2 , ■ ■ ■ , t n are the measurements of viral load 
and CD4+ T cell count from clinical studies which are subject to measure- 
ment errors. 

Let V(t), f(t), V'(t) and f'{t) be the estimates of V(t), T(t), V'(t) and 
T'{t) from Section 2.1, respectively, and substitute these estimates in (2.1). 
We then obtain a linear regression model [Liang and Wu (2008)], 

(2.2) V'{t) = a + atiT(t) + a 2 f'{t) - cV(t) + e(t), 

where e(t) includes all substitution errors. Using the least squares regression 
technique, we obtain the estimates of ao, a\ and a 2 which are termed as 
the pseudo- least squares (PsLS) estimates in Liang and Wu (2008). Notice 
that, from the above derivation, we have the relationship A = —ao/a 2 and 
p = a±/a 2 . Thus, we can obtain the estimates of c, A and p from model (2.2). 

Write a = (ao, ai, a 2 , c) and A = diag(l, T, T", — V). Let Pi{K) = 
J^_ 1 z l K(z)dz for I = 0,1,2. Using the arguments similar to those in Liang 
and Wu (2008), we can establish asymptotics for the smoothing-based esti- 
mators a of a. The Delta-method can then be used for establishment of the 
asymptotic properties of the estimators for A and p. Notice that, when the 
mean of e(t) is zero, the least squares estimates of ao, ct\ and a 2 are unbi- 
ased. However, the substitution error e(t) in model (2.2) is not mean zero, 
instead its mean is in the order of h 2 , and variance is of the order (n/i) _1 
which goes to zero when n — > oo. Thus, e(t) is different from a standard 
measurement error with mean zero and constant variance. 
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2.3. Stage III: Semiparametric regression for estimation of both constant 
and time-varying parameters. Now we propose how to estimate the remain- 
ing constant parameters, TV and 5, and the time- varying parameter, rj(t). 
Recalling (1.2), we have 

(2.3) jTj(t) = r,(t){T(t) - Tj(t)}V(t) - STj(t), 

while from (1.3) we can obtain 

A combination of (2.3) and (2.4) deduces 

y (2) (i) + cV'(t) = r)(t)N6T(t)V{t) - r)(t){V'(t)V(t) + cV 2 (t)} 

(2.5) 

-5{V'(t) + cV(t)}. 

Let Z(t) = V^{t] + cV'(t), Ui(t) = -{V'(t)V(t) + cV 2 (t)}, U 2 (t) = -{V'(t) + 

cV(t)}, U-s(t) = T(t)V(t), in which all functions and parameters have been 
estimated from Stages I and II. Thus, from (2.5), we can formulate a semi- 
parametric time- varying coefficient regression model [Wu and Zhang (2006)]: 

(2.6) Z(t) = U^S + U 2 (t)r](t) + U 3 (t)N6rj(t) + e*(t), 

where £*{t) includes all substitution errors. Note that 5 and N5 are con- 
stant parameters, while r/(i) is an unknown time-varying parameter. To fit 
this model, we can approximate the time- varying parameter rj(t) by the local 
polynomial approach, a basis spline method or other nonparametric tech- 
niques [Wu and Zhang (2006)]. The basis spline approach is straightforward 
and simple if the function rj(t) does not fluctuate dramatically as in our case 
of HIV dynamic models. Thus, we select to use the B-spline approximation 
here, that is, 

s 

(2.7) ri(t)^ aj b jtk {t), 

where a,- are constant B-spline coefficients, and bj t k(t) the basis functions 
of order k. For more details about B-splines such as construction of higher 
order basis functions via recurrence relations, the reader is referred to de 
Boor (1978). Thus, model (2.6) can be approximately written as 



Z(t) = U x {t)5 + J2ihk(t)U2(t)}aj 



(2.8) 



+ Y,{hk(t)U 3 (t)}N5a 3 + e*(t). 



PARAMETER ESTIMATION IN NONLINEAR ODE MODELS 



9 



This becomes a standard linear regression model and we can easily estimate 
the unknown constant parameters 5, a,j and N8a,j (j = 1, 2, . . . , s) from which 
we can derive the estimates of 6, N and ij(t). Note that AIC, BIC or AICc 
can be used to determine the order of the B-splines in our numerical data 
analysis. See detailed discussion on this in the next section. 

Also notice that Wu and Ding (1999) suggested that fitting viral dynamic 
data to different models for different time periods may result in better esti- 
mates of viral dynamic parameters. For example, we may fit a viral dynamic 
model with a closed-form solution to the early (first week) viral load data 
to obtain a better estimate of parameters 5 and c as proposed by Perelson 
et al. (1996). These estimates can then be substituted into the regression 
models in Stages II and III to obtain the estimates of other parameters. We 
will adopt this alternative strategy in our real data analysis. 

In summary, we are able to use a multistage approach to derive the esti- 
mates of all dynamic parameters in an HIV dynamic model. This approach 
only involves parametric and nonparametric/semiparametric regressions and 
the implementation is straightforward. The numerical evaluations of ODEs 
are avoided and the initial values of the state variables are not required. 
However, there are two issues that should be addressed for this approach. 
First, in this study, we focused on models (1.1)-(1.3) and the three stages 
were therefore particularly proposed for the specific model. However, the 
MSSB approach itself is a general method for estimating parameters in ODE 
models; for general guidelines in use of the MSSB approach, the interested 
reader is referred to Liang and Wu (2008). Second, although the estimators 
based on the MSSB approach have some attractive asymptotic properties, 
some limitations exist for this approach. In particular, it needs to estimate 
the derivatives of state variables which are sensitive to measurement noise 
when the data are sparse. This may result in biased estimates of parameters 
and measurement errors. In order to improve the estimates, we propose the 
spline-enhanced least squares approach in the next section. 

3. The spline-enhanced nonlinear least squares (SNLS) approach. In 

this section we introduce a spline-enhanced nonlinear least squares (SNLS) 
method to refine parameter estimates for models (1.1)-(1.3). The basic idea 
is to approximate the time- varying parameter rj(t) by the B-spline approach. 
Then models (1.1)— (1.3) become the standard ODEs with constant param- 
eters only. We can use the standard NLS approach to estimate the constant 
dynamic parameters and B-spline coefficients by numerically evaluating the 
ODEs repeatedly. This method is computationally intensive. Li et al. (2002) 
used a similar idea to approximate the unknown time-varying parameter in 
a pharmacokinetic ODE model. 



10 



H. LIANG, H. MIAO AND H. WU 



3.1. Spline approximation and nonlinear least squares. Different types 
of splines can be constructed based on different basis functions, such as the 
well-known piecewise polynomial splines and basis splines. Note that, for an 
arbitrary spline function of a specific degree and smoothness over a given 
domain partition, a linear combination of basis splines of the same degree 
and smoothness over the same partition can always be found to represent 
this spline function [de Boor (1978)]. Therefore, B-splines can be employed 
to approximate the time-varying parameter n(t) in this study without loss 
of generality. 

Similar to the approximation (2.7) in the previous section, n(t) can be 
approximated by a B-spline of order k with s control points [de Boor (1978)], 
that is, rf(t) ~ a jbj,k(t)- In addition, it should be noted that once the 
number and positions of control points are determined, the number and 
positions of knots are also automatically determined at the knot average sites 
[de Boor (1978)]. Thus, our model equations (1.1)-(1.3) can be approximated 
by 

(3.1) jTu(t) = A - P Tu(t) - { J>A*(*) }Tu(t)V(t), 



(3.2) jTjit) = |£aA*(t)JzM*M*) ~Sm), 

(3.3) jV(t) = N5Tj(t)-cV(t), 

where = (A, p, N, S, c, oi, 02, . . . , a s ) T are unknown constant parameters. 
Note that we have measurements of total CD4+ T cell counts, T(tj) = 
Tjj(U) + Tj(ti), and viral load, V(ti), that is, the measurement model can 
be written as 

(3.4) Y li = T(t i )+e li , i = l,2,...,n T , 

(3.5) Y 2j = V(t j ) + e 2j , i = l,2,...,ny, 

where En and £2j are assumed to be mean zero with constant variances 
and independent. Then the standard NLS estimator can be obtained by 
minimizing the objective function 

tit Tiy 

(3.6) RSS(0) = J2{Y U - T(ti, 6)} 2 + £{F 2i - V(tj,6)} 2 , 

i=i j=i 

where ut is the total number of CD4+ T cell measurements and ny is the 
total number of viral load observations; and T(ti,0) and V(tj,6) are numer- 
ical solutions to equations (3.1)-(3.3) using the Runge-Kutta method. How- 
ever, if en and £2j are correlated, the weighted (generalized) NLS method 
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can be used. To stabilize the variance and the computational algorithms, 
the log-transformation of the data is usually used in practice. More gener- 
ally, the weighted NLS approach can be used to more efficiently estimate 
the unknown parameters if the weights for different terms in the objective 
function are known. Usually the Fisher-Information-Matrix can be used to 
obtain the confidence intervals for the unknown parameters, but the boot- 
strap approach is more precise although it is more computationally intensive 
[Joshi, Seidel-Morgenstern and Kremling (2006)]. We will use the bootstrap 
method in our real data analysis. 

3.2. Hybrid optimization and spline parameter selection. For the SNLS 
approach, a critical step is to minimize a multimodal objective function 
(3.6) over a high-dimensional parameter space. In practice, it is challenging 
to find the global solution to such problems if the parameter space is high- 
dimensional (as always the case in many biomedical problems), and/or the 
parameter values are of different orders of magnitude, and/or the objective 
function is multimodal or even not smooth with noisy data. Thus, it is 
critical to develop an efficient and stable optimization algorithm. 

There are three main categories of optimization methods: direct search 
methods, gradient methods and global optimization methods. However, both 
the direct search methods (e.g., the Simplex method) and the gradient meth- 
ods (e.g, the Levenb erg-Mar quardt method and the Gauss-Newton method) 
can be easily trapped by local minima and even just fail for our problems 
[Miao et al. (2008)]. For details and application of such methods in ODE 
models, the reader is referred to Nocedal and Wright (1999) and Engle- 
zos and Kalogerakis (2001). Therefore, the global optimization methods are 
more suitable for the parameter estimation problem for ODE models. Moles, 
Banga and Keller (2004) compared the performance and computational cost 
of seven global optimization methods, including the differential evolution 
method [Storn and Price (1997)]. Their results suggest that the differen- 
tial evolution method outperforms the other six methods with a reasonable 
computational cost. Unfortunately, existing global optimization methods are 
very computationally intensive. Improved performance can be achieved by 
combining gradient methods and global optimization methods, called hy- 
brid methods. A hybrid method based on the scatter search and sequential 
quadratic programming (SQP) has been proposed by Rodriguez-Fernandez, 
Egea and Banga (2006), who showed that the hybrid scatter search method 
is much faster than the differential evolution method for a simple HIV ODE 
model. However, our preliminary work did not show improved performance 
of the scatter search method without SQP with respect to the differen- 
tial evolution method in terms of computational cost and convergence rate. 
Thus, the performance improvement of the hybrid scatter search method is 
mainly due to the incorporation of the SQP local optimization method. In 
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addition, Vrugt and Robinson (2007) suggested that a significant improve- 
ment of efficiency can be achieved by combining multiple global optimiza- 
tion algorithms. Here we combined the differential evolution and the scatter 
search method and incorporated with the SQP local optimization method 
[e.g., SOLNP by Ye (1987)] for parameter estimation of ODE models. We 
present the details of the proposed algorithm in the Appendix. 

Another critical problem that needs to be addressed in order to fit mod- 
els (3.1)-(3.3) is to determine the order k and the number and positions of 
control points s for the spline approximation. In general, we assume that 
rj(t) is a first-order continuous function of time. Therefore, B-splines of order 
2 (piecewise straight lines) should not be considered. Also, since the high or- 
der B-spline approximation (e.g., k > 5) may introduce unnecessary violent 
local oscillations (called Runge's phenomenon) [Runge (1901)], we consider 
up to 4th order B-splines in this study. As to the positions of the control 
points, to account for highly-skewed data, we select the control points' posi- 
tions such that they are equally-spaced in the logarithm time scale. We can 
use AIC, BIC and AICc criteria [Akaike (1973), Schwarz (1978), Burnham 
and Anderson (2004)] to determine the order and the knots, that is, 

(3.7) AIC = -2lnL + 2K, 

(3.8) BIC = -2lnL + K]n(N), 

(3.9) jucc-AW + mZg. 

where L denotes the likelihood function, N the total number of observations 
and K the number of unknown parameters. Under the normality assumption 
of measurement errors, these model selection criteria can be rewritten as 

(3.10) AIC = Nln(^pj +2K, 

(3.11) BIC = Nln(^pj +Kln(N), 

where RSS is the residual of sum squares obtained from the NLS model 
fitting. It should also be noticed that, as the number of spline knots goes to 
large, the associated parameter space will be high-dimensional. In this case, 
the novel and efficient optimization methods such as the hybrid optimization 
algorithms are necessary to locate the global minima of the NLS objective 
function. In addition, it is necessary to predetermine a good and informa- 
tive search range for each of the unknown parameters in order to ease the 
computational burden of the proposed optimization algorithm. Thus, it is 
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necessary and a good strategy to combine the proposed two approaches to 
estimate the dynamic parameters in a complex dynamic system. That is, 
first the MSSB approach introduced in the previous section can be used to 
obtain a rough estimate and the search range for the unknown parameters, 
and then the SNLS approach can be used to refine the estimates. 

4. Simulation studies. In this section we evaluate the performance of the 
MSSB and the SNLS approaches based on equations (1.1)— (1.3) by Monte 
Carlo simulations. An ad hoc bandwidth selection procedure is used for se- 
lecting a proper bandwidth h for MSSB. See Liang and Wu (2008) for more 
detailed discussions. The following parameter values were used to generate 
simulation data: T v (0) = 600, T/(0) = 30, V(0) = 10 5 , A = 36, p = 0.108, 
N = 1000, 5 = 0.5, c = 3, and r/(i) = 9 x 10~ 5 x {1 - 0.9cos(vrt/1000)}. Let 
T = Tjj + Tj denote the total number of infected and uninfected CD4+ T 
cells and V denote the viral load, the measurement models (3.4)— (3.5) were 
used to simulate the observation data with measurement noise, where £u and 
E2j were assumed to be independent and generated from normal distribu- 
tions with mean zero and constant variances o\ and o\, respectively. Equa- 
tions (1.1)— (1.3) were numerically solved within the time range [0, 20] (days) 
using the Runge-Kutta method, and solutions were output for equally- 
spaced time intervals of 0.1 and 0.2 which correspond to the number of mea- 
surements 200 and 100, respectively. The measurement errors a\ = 20,30,40 
and o\ = 100, 150, 200 were added to the numerical results of the ODE model 
according to the measurement models (3.4)-(3.5), respectively. 

To evaluate the performance of the MSSB and SNLS approaches for 
smaller samples sizes and larger variances that are similar to the actual 
clinical data evaluated in the next section, we also performed the simu- 
lations for the number of measurements n = = ny = 30 and 50, and 
o\ = 20 2 ,30 2 ,40 2 and o\ = 50 2 , 75 2 , 100 2 . To evaluate the performance of 
the estimation methods, we use the average relative estimation error (ARE) 
which is defined as 



where 6j is the estimate of parameter from the jth simulation data set 
and N is the total number of simulation runs. We applied the MSSB and 
SNLS approaches to 500 simulated data sets for each simulation scenario. 
We used the AICc and determined that the time-varying parameter rj(t) was 
approximated by a spline of order 2 with 3 knots in our simulation studies. 

The simulation results are reported in Table 1. From these results, we can 
see that the SNLS approach significantly improves the performance of the 
MSSB approach in all simulation cases as we expected. When the number 
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Table 1 

Simulation study: a comparison of the MSSB and SNLS approaches. The average relative 
error (ARE) was calculated based on 500 simulation runs. The time-varying 
parameter ri(t) is approximated by a 3-knots spline of order 2 



ARE (%): MSSB ARE (%): SNLS 



n 






A 


P 


N 


S 


c 


A 


P 


N 


S 


c 


30 


400 


2500 


183.70 


277. .7 '4 


104.55 


110.60 


92.37 


14.50 


28.00 


5.95 


6.31 


1.16 




900 


5625 


345.03 


537.24 


105.48 


130.67 


93.26 


18.90 


36.40 


7.60 


8.38 


1.72 




1600 


10,000 


411.94 


668.88 


112.47 


139.69 


94.96 


22.80 


44.60 


9.34 


10.70 


2.19 


50 


400 


2500 


77.41 


99.53 


81.93 


308.40 


65.46 


11.70 


22.80 


4.66 


4.69 


0.68 




900 


5625 


98.32 


106.77 


104.55 


384.52 


66.59 


15.20 


29.40 


6.10 


6.41 


0.95 




1600 


10,000 


267.91 


259.13 


112.00 


458.14 


67.86 


19.30 


38.00 


7.68 


8.20 


1.18 


100 


20 


100 


24.37 


40.89 


15.70 


65.18 


36.92 


2.50 


5.89 


1.54 


1.16 


0.12 




30 


150 


25.13 


42.50 


16.37 


97.53 


37.65 


2.84 


6.44 


1.68 


1.14 


0.13 




40 


200 


27.56 


46.98 


16.45 


114.09 


38.84 


2.96 


7.19 


1.72 


1.18 


0.15 


200 


20 


100 


8.01 


14.17 


9.72 


153.74 


15.15 


1.59 


4.55 


1.52 


0.61 


0.09 




30 


150 


8.24 


14.52 


10.42 


156.17 


15.54 


2.32 


4.97 


1.06 


0.84 


0.10 




40 


200 


8.85 


15.62 


10.76 


156.83 


16.25 


2.13 


5.37 


1.32 


0.96 


0.12 



of measurements is large (e.g., n = ut = ny = 100 or 200) and variances are 
small, the SNLS approach performs extremely well, but the AREs of the 
MSSB estimates of some parameters such as p and S are very large. When 
the number of measurements is small (similar to our real data analysis in the 
next section, that is, n = n? = ny = 30 or 50) and the variances are large, the 
performance of the SNLS estimates is still reasonably good, but the MSSB 
approach performs poorly for all parameter estimates. Similarly, the estimate 
of time- varying parameter r](t) from the SNLS approach outperforms that of 
the MSSB approach. Thus, the simulation results and our experience suggest 
that it is a good strategy to use the MSSB approach to obtain rough search 
ranges for unknown parameters, and then use the SNLS approach to refine 
the estimates. 

5. Viral dynamic parameter estimation for AIDS patients. We applied 
the proposed MSSB and SNLS approaches to estimate viral dynamic param- 
eters from data for individual AIDS patients based on model equations (1.1)- 
(1.3). Two HIV-1 infected patients were treated with a four-drug antiretro- 
viral regimen in combination with an immune-based therapy. Frequent viral 
load measurements were scheduled at baseline and after initiating the com- 
bination treatment: 13 measurements during the first day, 14 measurements 
from day 2 to week 2, and then one measurement at each of the following 
weeks, 4, 8, 12, 14, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 64, 74 and 76, 
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respectively. Total CD4+ T cell counts were monitored at weeks 2, 4 and 
monthly thereafter. 

As suggested in Section 2, first we applied the nonlinear regression model 
in Perelson et al. (1996) to fit the viral load data for the first week to 
obtain the estimates of 5 and c. The estimation results for the two patients 
are as follows: Patient I, 8 = 1.09 and c = 2.46; and Patient II, 5 = 0.43 
and c = 3.78. Then we applied the proposed MSSB and SNLS approaches 
to estimate all other viral dynamic parameters including the time-varying 
parameter rj(t) in models (1.1)— (1.3). The log 10 transformation of the data 
was used in order to stabilize the variance and computational algorithms. 

A rough estimate and a reasonable search range for each of the unknown 
parameters and the initial values of the state variables were obtained using 
the MSSB method. Then we applied the SNLS approach to refine the esti- 
mates. For the SNLS approach, we selected the smoothing parameters using 
the AIC, BIC and AICc criteria. We considered all the combinations of two 
spline orders (3 and 4) and 8 different numbers of control points (from 3 
to 10), and the results are reported in Tables 2 and 3 for Patients I and II, 
respectively. Note that as a practical guideline, if the number of unknown 
parameters exceeds n/40 (where n is the number of measurements), the 
AICc instead of AIC should be used. For our clinical data, n is about 40, 
and the number of unknown parameters varies between 8 and 15 for differ- 
ent models (the unknown initial conditions were also considered as unknown 

Table 2 

Model selection for Patient I: the time-varying parameter r](t) is approximated 
using splines of order 3 or 4 with the number of knots from 3 to 10 



Model 


Spline order 


Number of knots 


AIC 


BIC 


AICc 


1 


3 


3 


-165.7 


-151.5 


-161.7 


2 




4 


-163.8 


-147.7 


-158.8 


3 




5 


-174.4 


-156.5 


168.4 


4 




6 


-173.6 


-154.0 


-165.6 


5 




7 


-172.9 


-151.4 


-162.9 


6 




8 


-177.0 


-153.8 


-165.0 


7 




9 


-177.1 


-152.2 


-163.1 


8 




10 


-173.9 


-147.1 


-156.9 


9 


4 


3 








10 




4 


-170.5 


-154.4 


-165.5 


11 




5 


-173.7 


-155.8 


-167.7 


12 




6 


-173.4 


-153.7 


-165.4 


13 




7 


-171.5 


-150.1 


-161.5 


14 




8 


-172.3 


-149.1 


-160.3 


15 




9 


-175.5 


-150.6 


-161.5 


16 




10 


-174.0 


-147.3 


-157.0 



16 H. LIANG, H. MIAO AND H. WU 

Table 3 

Model selection for Patient II: the time-varying parameter rj(t) is approximated 
using splines of order 3 or 4 with the number of knots from 3 to 10 



Model Spline Order 


Number of knots 








1 3 


3 


-246.5 


-229.1 


-244.5 


2 


4 


-244.5 


-225.0 


-241.5 


3 


5 


-250.7 


-229.0 


-246.7 


4 


6 


-245.0 


-221.1 


-241.0 


5 


7 


-252.4 


-226.3 


-246.4 


6 


8 


-247.3 


-219.1 


-240.3 


7 


9 


-249.5 


-219.0 


-241.5 


8 


10 


-247.9 


-215.3 


-238.9 


9 4 


3 








10 


4 


-245.7 


-226.2 


-242.7 


11 


5 


-248.5 


-226.7 


-244.5 


12 


6 


-243.9 


-220.0 


-239.9 


13 


7 


-249.8 


-223.7 


-243.8 


14 


8 


-245.7 


-217.4 


-238.7 


15 


9 


-247.1 


-216.6 


-239.1 


16 


10 


-244.0 


-211.4 


-235.0 



parameters), which is much larger than n/40 = 40/40 = 1. So the AICc is 
more appropriate for our case. In fact, the AICc converges to AIC as the 
number of measurements gets larger, thus, the AICc is often suggested to 
be employed regardless of the number of measurements [Burnham and An- 
derson (2004)]. Therefore, our model selection is mainly based on AICc, 
although the AIC and BIC scores are also reported in Tables 2 and 3 for 
comparisons. Prom Table 2, we can see that both AICc and BIC selected 
the best spline model for rj(t) as order 3 and 5 control points for Patient 
I, although the AIC selected a different model. From Table 3, we can see 
that the AICc selected the same spline model (order 3 with 5 control points) 
for Patient II and both AIC and BIC ranked this model as the second best 
model. Based on above discussions, we are confident that the spline model 
with order 3 and 5 control points is the best approximation to the time- 
varying parameter rj(t). Thus, we will mainly report and discuss the results 
from this model. 

Figure 1 shows the data and model fitting results from the best SNLS 
approach for the two patients. The fitted curves are reasonably well. Table 4 
reports the estimation results from both MSSB and SNLS methods for the 
two patients. The MSSB approach provided good search ranges of unknown 
parameters for the SNLS approach. From the SNLS estimates, the prolif- 
eration rates of uninfected CD4+ T cells are A = 397 and 45 cells per day, 
respectively for Patients I and II; the death rates are p = 0.49 and 0.10 per 
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(a) Patient I: fitting CD4+ T cell count (b) Patient I: fitting viral load data 
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Fig. 1. The SNLS model fitting for two AIDS patients. 
Table 4 

The estimated values and associated 95% confidence intervals ( CI) of viral dynamic 
parameters for two AIDS patients. Perelson's model [Perelson et al. (1996)] was first 
fitted to the viral load data within the first week to obtain estimates of 5 and c. For 
Patient 1,8= 1.09 and c = 2.46; for Patient II, 5 = 0.43 and c = 3.78. The time varying 
parameter rj(t) of both Patients I and II was approximated using a spline of order 3 with 

5 control points 



Parameter 



Patient I 



Patient II 



MSSB 
(CI) 



SNLS 
(CI) 



MSSB 
(CI) 



SNLS 
(CI) 



A (cell per day) 254.49 397.09 18.38 45.45 

(181.11, 327.87) (216.43, 594.30) (3.50, 33.27) (29.78, 81.48) 

p (per day) 0.34 0.49 0.04 0.10 

(0.24, 0.44) (0.26, 0.75) (0.00, 0.08) (0.06, 0.18) 

N (virion per cell) 1178.23 264.74 2769.32 1114.37 

(558.20, 1798.26) (203.40, 350.00) (2484.54, 3074.10) (856.62, 1428.93) 
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day which correspond to the half-life of 1.4 and 6.9 days for the two patients, 
respectively; the numbers of virions produced by each of the infected cells 
are 265 and 1114 per cell for Patients I and II, respectively. 

Note that the nonlinear regression estimates of the death rate of infected 
cells, 6, for the two patients are 1.09 and 0.43 per day with the corresponding 
half-life of 0.64 and 1.61 days, respectively, which indicates a higher death 
rate of infected cells compared to the uninfected cells for both patients. The 
estimates of viral clearance rate (c) were 2.46 and 3.78 per day with the 
corresponding half-life of 0.28 and 0.18 days for the two patients, respec- 
tively. These results are similar to those from the previous studies [Perelson 
et al. (1996, 1997)] and biologically reasonable. However, to our knowledge, 
this is the first time the estimates of the proliferation rate and death rate 
of uninfected CD4+ T cells (A and p) and the number of virions produced 
by an infected cell (iV) have been obtained directly from the clinical data 
of AIDS patients. From the estimation results in Table 4, we can see that 
the death rate of infected cells is 2 to 4 fold higher than that of uninfected 
cells for Patients I and II, respectively. We can also see that the difference in 
the estimated dynamic parameters between the two patients is large. This 
is consistent with the argument that the between-subject variation of viral 
dynamics in AIDS patients is large [Wu and Ding (1999), Wu et al. (1999)]. 
Thus, the personalized treatment is necessary for AIDS patients. 

The estimated time- varying parameter, the infection rate rj(t) for Patients 
I and II are plotted in Figure 2. From the estimates, we can see that the 
infection rates for the two patients have similar patterns which are primarily 
due to similar patterns of viral load and CD4+ T cell counts for the two 
patients. It seems that there was an initial fluctuation in the infection rate 
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Fig. 2. The estimated infection rate r)(t) (solid line) and its bootstrap 95% confidence 
interval (dashed line). 
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at the beginning of treatment, and then the infection rate was stabilized 
after one or two weeks. We can also see that the estimated infection rate is 
not zero, which may suggest the imperfectness of the long-term treatment 
and the development of drug-resistant mutants. This is inconsistent with 
the perfect treatment assumption in Perelson et al. (1996, 1997), although 
it may be valid in a short period of time after initiating the treatment. 

6. Discussion and conclusion. We have proposed two approaches to iden- 
tify all dynamic parameters including both constant and time-varying pa- 
rameters in an HIV viral dynamic model which is characterized by a set of 
nonlinear differential equations. This is a very challenging problem in the 
history of HIV dynamic studies. The proposed multistage smoothing-based 
(MSSB) approach is straightforward and easy to implement since it does 
not require numerically solving the ODEs and the initial values of the state 
variables are not needed. However, one limitation of this approach is that 
it requires the estimates of derivatives of the state variables, which are usu- 
ally poor when the data are sparse. This may result in poor estimates of 
unknown dynamic parameters. The second method that we proposed is the 
spline-enhanced nonlinear least squares (SNLS) approach. This method is 
more accurate in estimating the unknown parameters, but it requires nu- 
merical evaluations of ODEs repeatedly in the optimization procedure and 
the convergence of the computational algorithm is problematic when the pa- 
rameter space is high-dimensional. We have proposed a hybrid optimization 
technique to deal with the nonconvergence and the local solution problems, 
but the computational cost is high and a good search range for each of the 
unknown parameters is required. In practical implementation, we propose 
to combine the two approaches, that is, first using the MSSB approach to 
obtain rough estimates and possible ranges for the unknown parameters 
and then using the SNLS approach to refine the estimates. We have applied 
this strategy to estimate all dynamic parameters from data for two AIDS 
patients. 

It is very important to estimate all dynamic parameters in HIV dynamic 
models directly from clinical data when prior knowledge is not available or 
fixing some parameters will significantly affect the estimates of other pa- 
rameters in a dynamic model. To our knowledge, this is the first attempt to 
propose a procedure to estimate both constant and time- varying parameters 
in the proposed HIV dynamic model directly from clinical data. Although 
Huang, Liu and Wu (2006) have attempted to do so using a Bayesian ap- 
proach, strong priors were used for most of the dynamic parameters and 
a parametric form for the time-varying parameter was employed. Notice 
that the time-varying parameter rj(t) is a function of antiviral treatment 
effect in the HIV dynamic model [Huang, Liu and Wu (2006)]. It is very 
important to estimate this time-varying parameter for each AIDS patient 
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individually in order to better design the treatment strategy and person- 
alize the treatment for each individual. In addition, to better assess the 
treatment efficacy, the time-varying infection rate r/(t) can be further mod- 
eled as rj(t) = r] pre x (1 — r]treat{t)), where rj pre denotes the infection rate at 
the pre-treatment equilibrium, and r]t re at(t) the time- varying infection rate 
after treatment. By comparing r] pre with rj(t), it will be easy to quantify the 
treatment effects on infection rate. However, since r] pre is strongly correlated 
with rjtreat{t), rj pre has to be fixed to estimate r] tre at{t)- Unfortunately, the 
pre-treatment equilibrium data were not collected in this study such that 
we could not determine the value of rj pre and therefore to do the compari- 
son. Finally, the treatment could also affect parameter N; however, limited 
by the clinical measurements, we simplified our model without considering 
parameter N as time-varying. 

We used both viral load and CD4+ T cell count data in order to identify 
all dynamic parameters in the HIV dynamic models (1.1)— (1.3). Usually 
AIDS investigators believe that CD4+ T cell count data are very noisy with 
a large variation due to both measurement errors and natural variations 
of CD4+ T cell counts. Thus, it is important to improve the data quality 
in order to get more reliable estimates of dynamic parameters, although it 
is beyond statisticians' control. Our study also suggests that the frequent 
measurements, in particular, the viral load measurements during the early 
stage after initiating an antiviral therapy, are important to identify some 
constant dynamic parameters such as c and 5, and long-term monitoring of 
viral load and CD4+ T cell counts is necessary to estimate other parameters, 
in particular, the time-varying parameter. 

APPENDIX: DIFFERENTIAL EVOLUTION AND SCATTER SEARCH 

The differential evolution (DE) algorithm [Storn and Price (1997)] is a 
typical evolutionary algorithm that searches the optimum by inheritance, 
mutation, selection and crossover of parent populations (e.g., vectors of pos- 
sible parameter values in the parameter space). The initial population is 
randomly generated from a uniform distribution within the search range to 
cover the whole region, denoted by Xi t c {i = 1, 2, . . . , Np), where Np is the 
number of members in this generation. The subsequent population inher- 
its and mutates by randomly mixing the previous generation with certain 
weights. Storn and Price (1997) recommended a mutation method 

(A.l) Vi,G+l = x n ,G + F(x r2)G - x r3) a), 

where v^G+i is the ith member of generation (G + 1), x ri) a, ^r 2 ,G and x T3i g 
the members of the previous generation G, and F > the amplification 
factor. Integers r\ , T2 and are randomly chosen from {1,2,..., Np} which 
are mutually different from each other and different from i. Also, to increase 
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diversity, the mutated member i^g+i exchanges its components with Xi t G 
with a given probability, the crossover ratio. For this purpose, a number 
within [0, 1] is generated for each component of ^i,G+i by a uniform random 
number generator. If this number is greater than the crossover ratio, then 
the component of Xi t a is kept. Finally, the best member in the mutated 
population is selected by comparing the values of the objective function. 
The convergence rate of the differential evolution method depends on specific 
mutation strategies used, that is, the amplification factor and the crossover 
ratio. This method has been successfully employed in previous studies [Miao 
et al. (2008, 2009)]. For more details about the DE algorithm, the reader is 
referred to Storn and Price (1997). 

The scatter search method was first proposed by Glover (1977) and ex- 
tended later by Laguna and Marti (2003, 2005). The scatter search method 
is not a genetic algorithm; instead, it locates the optimum by tracking the 
search history and balancing visiting frequency within different segments of 
the search region using sophisticated strategies. Here we give a brief intro- 
duction to this method. For more details, the reader is referred to Laguna 
and Marti (2003) and Rodriguez-Fernandez, Egea and Banga (2006). 

Let li and Ui denote the lower and upper boundary of the ith component 
of the parameter vector, respectively. Then the region between Zj and Ui can 
be divided into to segments (e.g., m = 4), which are usually of equal length. 
Let Sij (j = 1, 2, . . . , to) denote the jth segment of the search region [li,Ui] . 
If a parameter value candidate of 0i falls in s^j , we call it a visit to s^j . Let 
fij denote the total number of visits to Sy , then a visiting history matrix F 
can be formed: 



Since the scatter search method is population-based, the first population 
needs to be generated to start the algorithm. Let Nq denote the total number 
of parameter vectors in the first population, then to parameter vectors are 
generated first such that the ith component of the jth vector falls into Sj,-. 
In this way, all elements of matrix F become one. To generate the rest 
(Nq — to) parameter vectors in the first population, a random vector z = 
(zi, Z2, ■ ■ • , z q ) T is generated for each parameter vector with Zi (i = 1,2, ... ,q) 



(A.2) 



F = 



fll fl2 

hi Hi 





k=l 1 /Jik 
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following a uniform distribution on [0,1]. Then the zth component of the 
parameter vector to be generated falls into Sit for the smallest k satisfying 

k 

Zi<^Pij, k = l,2,...,m. 
i=i 

Note that after each new parameter vector is generated, the matrix F is also 
updated. 

Once the first population is generated, the parameter vectors are divided 
into two categories: elite vectors and diverse vectors. Let n e denote the 
number of all elite vectors, which is usually a small fixed number (e.g., 20). 
Then n e /2 elite vectors are chosen based on the smallest objective function 
values; the rest half of the elite vectors are chosen to be farthest from all the 
first half elite vectors in the sense of Euclidean distances. Thus, both fitness 
and diversity are considered in the construction of the elite vectors. Then the 
new parameter vectors can be generated by combination of the elite vectors. 
Let x^ l > and denote two different elite vectors and has a smaller 
objective function value than x^ 2 ) . Then three types of combinations can be 
employed to generate new parameter vectors: 

• Type 1: p\ = x^ — d, 

• Type 2: p2 = x^ 1 ' + d, 

• Type 3: p 3 = x^ + d, 

where d = r T ■ (x^ — x^) with all the components of r generated from an 
uniform distribution on [0, 1]. If both x^ 1 ) and x^ are in the first half elite 
vectors in terms of fitness, one new vector of type 1, one of type 3 and two 
of type 2 are generated; if only x^ belongs to the first half elite vectors 
in terms of fitness, one new vector of each type is generated; if both x^ 
and x^ 2 ) belong to the second half elite vectors in terms of farthest distance, 
then one new vector of type 2 and one of either type 1 or 3 are generated. 
The fitness of these new generated vectors is then compared to that of the 
elite vectors. The new vectors with smaller objective function values (that 
is, better fit) will replace the elite vectors with larger objective function 
values. This procedure is repeated until no new vectors can replace any 
elite vectors. Now the algorithm can either stop or continue by regenerating 
diverse vectors. 
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